Modelling of quantum information processing with Ehrenfest guided trajectories: 

a case study 
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We apply a numerical method based on multi-configurational Ehrenfest trajectories, and demon- 
strate converged results for the Choi fidelity of an entangling quantum gate between two two-level 
systems interacting through a set of bosonic modes. We consider both spin-boson and rotating wave 
Hamiltonians, for various numbers of mediating modes (from 1 to 100), and extend our treatment 
to include finite temperatures. Our results apply to two-level impurities interacting with the same 
band of a photonic crystal, or to two distant ions interacting with the same set of motional degrees 
of freedom. 
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I. INTRODUCTION 

The ability of tracking the evolution of complex quan- 
tum systems will be a crucial support to the design 
and development of future quantum technologies. A 
paradigm of particular interest for the latter is one where 
finite dimensional quantum systems, typically two-level 
systems (qubits), interact through a 'bus', made up of 
a set of bosonic field modes (ll— [331] . In this paper, we 
will consider the case of two qubits interacting with a 
common discrete set of modes in the non-perturbative 
regime. While the possibility of creating entanglement, 
even at steady state, by the interaction with a common 
bosonic bath has been highlighted repeatedly in the lit- 
erature, the case of non-perturbative interactions with 
a bunch of 10 or 20 modes (as will be the case in our 
study) presents several major technical difficulties, es- 
sentially due to the impossibility of an analytic master 
equation approach - only possible in the continuum limit 
under the Born-Markov approximation - and to the huge 
size of the dynamically relevant part of the Hilbert space. 

Various numerical approaches have been developed to 
emulate these dynamics on classical computers, such 
as the so called multiconfigurational time dependent 
Hartree method [3^ - l36j and its "Gaussian" variation 
37 1. various schemes based on path integral techniques 
38ll39|. and even the adaptation of time-adaptive density 
matrix renormalisation group techniques [40| borrowed 
from condensed matter theory. Here, we will tackle such 
difficulties by borrowing a method co-pioneered and de- 
veloped by one of the authors in the arena of chemical 
physics [4l|, |42[ • The method is based on the adoption 
of a set of tensor products of time-dependent coherent 
states as a discrete 'basis grid' on which to represent the 
field degrees of freedom (referred to as "coupled coherent 
states" in the literature HH), and on letting the states 
of the co-moving grid evolve according to their Ehrenfest 
dynamics (whose application to a grid of coherent states 
goes under the name of "multi-configurational Ehrenfest" 
method). This approach has the advantage of being rel- 
atively light in terms of computational resources, easy 



to program, and yet of allowing one to follow coherent 
quantum dynamics in detail, as it will be shown. 

In our study, we will focus on a specific, but very rele- 
vant, aspect of the quantum dynamics of the two qubits: 
we shall consider the realisation of an entangling quan- 
tum gate between them, namely of a controlled Z (CZ) 
gate. To estimate the quality of such a realisation we 
will consider the quantum fidelity between the pure state 
corresponding to the CZ gate by the standard channel- 
state duality (Choi isomorphism (44l - |4r| ) and the quan- 
tum state corresponding to the channel acting on the 
two qubits upon partial tracing over the field's degrees of 
freedom. 

Our main aim is demonstrating the capability of 
Ehrenfest guided trajectories in phase space to produce 
reliable and converged results for complex figures of 
merit, able to reveal detailed information about the quan- 
tum dynamics of the constituents. The 'Choi fidelity' of 
a two-qubit quantum gate is a property of the dynam- 
ics itself, and not of the initial state, and its evaluation 
requires, at any time, the evolution of ten initial states: 
it is, therefore, a rather cumbersome figure of merit to 
compute, let alone to optimise over a rather large range 
of values of the dynamical parameters, as we did. The 
advantages of the Ehrenfest guided trajectories over - ar- 
guably more precise but heavier - approaches based on 
full variational principles are manifest in such circum- 
stances. 

As for direct impact, let us remark that our study 
would apply on systems like two-level impurities trapped 
in a photonics crystal and interacting with the same band 
of allowed modes [HI, |U EH , or to the internal levels of 
two ions interacting with all the vibrational modes of an 
array of ions in a linear trap [49[ . It is important to point 
out that our treatment can account for finite, although 
relatively small, temperatures as well. 

Our paper is organised as follows. In section [TT| we will 
review the basic theory behind methods of solution of the 
Schrodinger equation based on a set of time-dependent, 
Ehrenfest guided basis states. We will not dwell so much 
on the technical details, which are covered elsewhere, as 
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on the basic concepts, and will try to present them in 
terms which will be friendly to an audience with no pre- 
vious familiarity with the terminology of chemical physics 
or molecular dynamics. In section IIIII we will introduce 
the physical Hamiltonian and define precisely our chosen 
figure of merit. Section ITVl will contain the main results 
of our numerical study. Finally, we will draw conclusions, 
and discuss advantages and shortcomings of our method 
of choice, in section IVTl 



II. EHRENFEST GUIDED TRAJECTORIES 

The main difficulty in dealing with a system including 
few two-level systems and a bunch of M bosonic modes 
is clearly how to handle the infinite dimensional bosonic 
Hilbcrt space. The method we apply here, referred to 
in the literature as 'multi-configurational Ehrenfest' and 
first introduced in 4l| , tackles this difficulty on the shoul- 
ders of two major assumptions: 

i) the state space of the field modes is represented as a 
superposition of N time dependent coherent states; 

ii) the time dependence of the coherent states is deter- 
mined by a simplified variational principle (which, 
in other words, dictates how the finite dimensional 
subspace spanned by the set of coherent states 
changes with time, trying to keep it in the dynam- 
ically relevant region). 

The finite dimensional systems involved are, on the other 
hand, treated by representing their entire Hilbert space, 
spanned by d orthonormal basis states 1 I), for I € {1, d}. 
In agreement with i), the ansatz wave- function of the 
whole system reads 



JV 



(i) 



i=i j=i 



coherent states: | otj) = ®„ = i | &j) 



where otj £ <C V j and each | otj) is a tensor product of 

such that, if a m is 
the annihilation operator of mode m, one has a m \ otj) = 

aj™''| OLj). Since we will be dealing with two qubits, it 
will be d = 4 for us. 

The evolution of the dynamical parameters is more 
conveniently described by adopting a Lagrangian formu- 
lation. For a Hamiltonian operator H , let us define C 
as 



C = {i)\H -id t \i>) 



(2) 



The time-derivative operator is defined as the differentia- 
tion of the time dependents coefficients c; j and by the re- 
lationships dt\l) = (the finite dimensional system's ba- 
sis is time-independent) and d t \ otj) = Y^m=l\-®j ( a m ~ 
a-™^* /2) — /2]\a.j) (derived from the time- 

dependence of a coherent state with varying phase space 



position) . The quantity C is hence a function of the co- 
efficients c/j, the complex parameters ctj and their time- 
derivatives bi j and a.j. In fact, it can be shown that C 
serves as a Lagrangian for the quantum system, in the 
sense that the Euler-Langrange equations 



dC 



d dC 



dc. 



dt dc 



1,3 



(3) 



are equivalent to Schrodingcr equation [68j. See Ap- 
pendix [B] for more details. 

Besides determining the state evolution, the varia- 
tional formalism also provides one with a recipe to update 
the basis such that the expression C of Eq. ((2]), which 
clearly always equals in the exact dynamics, is min- 
imised during the time-evolution. Such a minimisation, 
which in essence keeps the basis in the 'most relevant' 
region of the Hilbert space within the constraints of the 
adopted approximation, would be obtained by consid- 
ering the full Euler-Lagrange equations for the M x N 
complex parameters ctj and their time-derivatives. This 
would be a large nonlinear system of coupled equations, 
requiring a substantial numerical effort to be solved. In- 
stead, we introduce here assumption ii), and replace the 
full variational equations for otj with a simplified ver- 
sion thereof. In particular, we will neglect all terms cou- 
pling the different aj's, on the grounds that the overlaps 
(ctj \otk) are typically very small if the number of modes 
M is large enough. For each j, let us then define the 
vector | ipj) as 



®\OLj 



(4) 



i=i 



and the corresponding 'approximated' Lagrangian Cj as 

C j = (i> j \H-id t \4, j ) . (5) 

The equation of motion for the parameter otj (the m-th 
component of the vector otj) is 



c)c\ 



d dCj 

dt gArr 



(G) 



where we also neglect the time-dependence of the coef- 
ficients c/j, such that each Lagrangian Cj only depends 
on the four complex parameters and on the M com- 
plex parameters represented by the entries of oij (and 
on their time-derivatives ctj, see Appendix [A] for fur- 
ther details). Eq. (JB]) defines the "Multi-Configurational 
Ehrenfest" (MCE) method we are using. 

Notice that the assumption ii) is not, per se, an ap- 
proximation, but rather just a way of choosing the time- 
dependence of the adopted basis. However, it should be 
stressed that, in general, the exact Euler Lagrange equa- 
tion for the full variation of the parameters otj is likely 
to provide one with a more accurate result in that it will 
yield a smaller Lagrangian C (which is zero in the exact 
dynamics). 
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However, Eq. ([6]) is much easier to treat numerically, 
hence the advantage of our method, which can be eas- 
ily programmed and applied with modest computational 
resources and often provides results in very good agree- 
ment with complete variational methods, like MCTDH 
or G-MCTDH jaUHzJ. 

Multi-configurational Ehrcnfcst guided trajectories 
have been thoroughly tested for spin-boson dynamics un- 
der different spectral densities, establishing the reliability 
of their converged results in several, diverse situations 
fill l5(il ]. Here, we will instead apply them to study a 
composite system including discrete sets of bosonic field 
modes, where coherent quantum information processing 
can be carried out. 



III. MODEL AND FIGURE OF MERIT 

We set out to study coherent quantum information 
processing for a system comprising two two-level sys- 
tems ( "qubits" ) connected by M bosonic modes through 
a spin-boson like coupling. 

In principle, this represents the archetype of a quan- 
tum system where complex dynamics and information 
processing tasks can be carried out, and whose dynamics 
is impervious to non-approximated methods. In prac- 
tice, our case study may be thought of as representing 
two two-level atoms (or impurities) interacting with the 
same photonic band of a photonic crystal (47| . or a sim- 
ulation of the same setting in a linear array of trapped 
ions (where the qubits are embodied by internal levels 
of the ions interacting with the same set of vibrational 
normal modes |49f). 

For future convenience, let us re-label the four states 
of the computational basis of the two two-level systems 
as follows: 



The operators &x , o'x 1 , a)?* and af, will stand for 
the customary Pauli operators in the Hilbert spaces of 
qubit 1 and 2. For instance, in the adopted basis, a y x 



1) 

2) 
3) 
4) 

J2) 



It) , 

U) , 

tt) . 



(7) 
(8) 
(9) 
(10) 



~(i) 



,(2) 



(1) 



is defined by a4 | 1 
.Ml 



3),ai 1} |3) 



1) 



^12) 



4), 



<ri'|4) = |2). 

In our study, we shall consider both an actual 'spin- 
boson like' Hamiltonian: 



H 



E2 ^ M 
7 — 1 



tt) 



+ 



tt) tt) 1 

1m u x X^r 



(ii) 



and its rotating wave counterpart: 



so 

+ 



tt) 



(12) 



where a 



tt) 



As well known, the 



Hamiltonian H rw is a good approximation of H in the 
almost resonant, high frequency case that is, in our nota- 
tion, for |2e — uj m \ <C |2e + w m |, Vm. In the following, we 
will consider systems with different numbers of bosonic 
modes M, various values of frequencies {w m } and spin- 
boson couplings {gm^}, and different parameters Ai and 
A-2- Also, we will set h = 1 throughout the paper. 

In reproducing the dynamics of the two qubits by treat- 
ing the field through multi-configurational Ehrenfest tra- 
jectories, we will aim at obtaining converged results for a 
figure of merit of interest in the study of quantum infor- 
mation processing, namely the fidelity with which an en- 
tangling controlled Z (CZ) gate can be realised for the two 
qubits through the mediating bosonic modes. In terms of 
the basis states of Eqs. (f7l [T0|) . a CZ gate is represented 
as a unitary Ucz leaving all the basis states invariant 
except for | 4) , which becomes — | 4) , that is 



Ucz\j) = f(j)\j) for 



1 < j < 4 



(13) 



where f(j) — 1 for j G {1, 2, 3} and f(j) = —1 for j = 4. 
The relevance of a CZ gate to quantum information pro- 
cessing stems from its being a maximally entangling gate 
which, combined with singlc-qubit unitarics, forms a uni- 
versal quantum set for gate based quantum computation 

a 

At zero temperature, the quantum operation T t we 
want to compare with the CZ unitary gate is defined as 
follows, in terms of a notional initial density matrix of 
the qubits g: 



T t (g) = Ti B 



-iHt 



(0® |0>(0|)e 



iHt 



(14) 



where Ttb stands for partial tracing over the Hilbert 
space of the bosonic modes and | 0) is the vacuum state 
of the modes. We will also extend our treatment to in- 
clude a finite temperature 1//3 of the bosonic modes (in 
natural units where ks = 1), in which case the quantum 
operation T t ,p will be given by 



T tt p(g) = Tt b 
with 



e~ lHt I g<E> 



f P p {a)\a)(a\d 2M a] 

(15) 



M 



771 — 1 



-l)|«m|- 



(16) 



The function Pp(at) is just the Glauber-Sudarshan P- 
representation of a thermal state of the bosonic modes 
(given by the product of individual P-representation for 
each of the modes). In our notation a <G (D 2M , while each 
component of a. is denoted by a m ■ Clearly, one has that 
lim^oo Yt.p = T t . 

In our numerical study, we will reproduce the oper- 
ations r t by adopting the method detailed in the pre- 
vious section, which is defined for pure states, and also 
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(a) 



(b) 
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FIG. 1: Choi fidelity F versus rescaled time, for H, w with 
e = A = 0, <?i = 1, g<i = 1.9, obtained at zero temperature 
by MCE method (dot-dashed) and exact analytic integration 
(dotted) for M = 1 and lo\ = 0.1 (a), and M = 3 and ui m — 
0.1m for 1 < m < 3 (b). The lines F = 0.25 are reported for 
reference. 



reconstruct the operations by sampling different ini- 
tial pure coherent states | a) for the field according to 
the distribution given by P(a). We will describe the 
field in terms of coupled coherent states during the time 
evolution and then trace it out to achieve the quantum 
operation acting on the two qubits. 

To define the gate fidelity F, we will make use of the 
classic channel-state duality (Choi isomorphism) map- 
ping linear quantum operations over a Hilbert space H 
into quantum states on the Hilbert space FL ® FL [45| . 
Turning to the two qubits Hilbert space FL spanned by 
the basis states (fT HTTJj) . let us define the maximally en- 
tangled fiducial state | ip) (belonging to FL 2 ) as: 

|tf) = i(|l>®|l> + |2)®|2> + |3)®|3) + |4)®|4». (17) 

For a generic CP-map f2, the corresponding quantum 
state gn may be defined as 

en = (fi®l)(|VW|), (18) 

where 1 is the identity map acting on FL. 

Since the CZ gate is unitary, the quantum state Qcz 
is bound to be pure: Qcz = | l Pcz)( l Pcz \ , with 

Wcz) = i(|l)®|l> + |2>®|2> + |3>®|3>-|4>®|4» 

= 5£/(i)(li>®li» • (19) 

The state gr t « corresponding to Tt,p is instead given by 

er*.e = \ E rt,p(\j){k\)®\j){k\. (20) 
i,k=i 

We can then naturally define the CZ operation fidelity F 
(which we shall informally refer to as 'Choi fidelity') as 
the overlap 

i 4 

F = (vcz \Qv t Jvcz) = Ye E 0' ^ (I ^ fc I) I fc ) ' 

j,k=l 

(21) 
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FIG. 2: MCE results for the Choi fidelity F versus rescaled 
time, for H rw with e = A = 0, gi = 1, M = 10, uj m = 0.1m 
for 1 < m < 10, zero temperature and different values of g2- 
The line F = 0.25 is reported for reference. 



The quantity F captures, in one real number, a rele- 
vant facet of the quantum dynamics governing the two 
qubits. Its relationship to quantum coherence is manifest 
in that if the off-diagonal elements between the basis vec- 
tors of Eqs. (jTl fTU)) are set to zero, then one has F < 1/4. 
Any value of F larger than 1/4 is thus in a sense a signa- 
ture of quantum coherence. More importantly, F is also 
a measure of how well a coherent quantum task can be 
performed and is also strictly related to the entanglement 
generated between the two qubits (in that entanglement 
a perfect CZ gate would get entanglement equal to 1 
ebit for a properly chosen initial state). Moreover, F, 
although partial to the chosen reference gate (CZ in this 
case) is completely independent of the initial state, and 
represents a property of the dynamics alone. 

Of course, we could have chosen more generic quanti- 
fiers like, for instance, the largest eigenvalue of the opera- 
tor gr t a j which would equal 1 in the ideal case where the 
qubits undergo a unitary evolution and would quantify, 
in a sense, the overall coherence of the qubits' evolution. 
However, we deem such choices to be less informative 
with regard to the applicative potential of a complex dy- 
namics. 

Given a potentially useful quantum dynamics, the 
knowledge of F is instead very desirable to possess. 
Demonstrating the use of a numerical technique capable 
of providing one with reliable estimates of F in relevant 
situations is, in a nutshell, the aim of the current analy- 
sis. 



IV. CHOI FIDELITY OF THE CZ GATE 

Here, we will slightly deviate from the previously 
adopted notation by setting = gj for all j and /. 
It is important to remark that assuming equal couplings 
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FIG. 3: MCE results for the Choi fidelity F versus rescaled 
time, for H rw with e = A = 0, g\ = 1, M = 10, u) m — 0.1m 
for 1 < m < 10, j3 = 10 and different values of 32- The line 
i 7 = 0.25 is reported for reference. 



FIG. 5: MCE results for the Choi fidelity F versus rescaled 
time, for H rw with e = A = 0, gi = 1, M = 10, uj m = 0.1m 
for 1 < m < 10, /3 — 5 and different values of g%. The line 
F = 0.25 is reported for reference. 
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FIG. 4: MCE results for the Choi fidelity F versus rescaled 
time, for H Tm with e = A = 0, gi = 1, M = 10, oj m = 0.1m 
for 1 < m < 10, /3 = 10 and different values of 32 (red stands 
for higher values, blue for lower values). 




FIG. 6: MCE results for the Choi fidelity F versus rescaled 
time, for H rw with e = A = 0, g\ = 1, M = 10, u) m = 0.1m 
for 1 < m < 10, P = 5 and different values of gi (red stands 
for higher values, blue for lower values). 



between each qubit and all the modes is in no way essen- 
tial to our numerical approach. Such an assumption can 
be - and will be, in the following - relaxed if needs be. 

Let us then, to begin with, set the coupling between 
the first qubit and the field modes g\ to 1, and essentially 
choose it as the unit of time. Let us also, until further 
notice, set e = A = 0, j3 — > 00 (zero temperature) and 
consider the rotating wave Hamiltonian H rw . Note that 
the Hamiltonian H rw with e — A — can be derived 
from the full Hamiltonian H with A = by switching to 
interaction picture and applying the rotating wave ap- 
proximation: e can be thus be set to zero and each field 
frequency ui m is shifted as per uj m — s- uj m — 2e. 

By exploiting the conservation of the number of excita- 
tions, the dynamics governed by the Hamiltonian H rw for 



A = can be easily solved analytically. The agreement 
between such analytical solutions and the MCE results 
has been tested for up to ten modes and is excellent. 



Figs. 1(a) and 1(b) show such an agreement in terms of 
CZ Choi fidelity F for 92 = 1.9 and, respectively, one 
mode with uj\ = 0.1 and three modes with uj\ = 0.1, 
L02 = 0.2 and W3 = 0.3. An initial peak with fidelity 
larger than 0.9 is immediately apparent: this peak will 
be the main object of our investigation, for larger num- 
bers of modes too. For three modes, the peak appears 
at a time which is approximately reduced by a y/3 factor 
with respect to the single mode case. This cooperative 
effect is confirmed for all number of modes up to 20, and 
is simply due to the fact that the qubits are coupled to 
the field through the mode -7^ X^m=i a m 5 w ith an ef- 
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FIG. 7: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H rw with e = A = 1, gi = 1, 
M = 10, uim = 0.1m for 1 < m < 10 and different values of 
(?2. The line F — 0.25 is reported for reference. 



FIG. 8: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H rw with e = A = 1, gi = 1, 
M = 10, Ldm = 0.1m for 1 < m < 10 and different values of 
g2 (red stands for higher values, blue for lower values). 



fective coupling which scales like vM (clearly, this is the 
consequence of assuming equal couplings with all modes) . 

Figs. EUS depict a detailed analysis of the Choi fidelity 
for M = 10 bosonic modes with frequencies uj m = 0.1m 
for 1 < m < 10, three different temperatures (/3 — > oo, 
/3 = 10 and (5 = 5), and different values of the coupling 
(?2, scanned over the range 1.8 — 2.7. The zero tempera- 
ture case (Fig. [2J shows how the dispersion of quantum 
coherence among the field's degrees of freedom affects 
the gate's fidelity (whose maximum is smaller than in 
the one- and three- modes cases), although in the consid- 
ered region of dynamical parameters the effect is not as 
pronounced as one could imagine. The plots clearly show 
the detrimental effect of thermal fluctuations on coherent 
quantum effects, even at such relatively small tempera- 
tures. In practice, this suggests severe constraints on the 
temperature for the observation of coherent quantum ef- 
fects mediated by discrete vibrational modes (typically 
much more susceptible to thermal excitations than opti- 
cal modes due to their lower frequency) , considering that 
the highest temperature accounted for is around 0.2<7i in 
natural units. Most importantly, we were able to deter- 
mine the optimal value of the coupling gi with respect to 
a vast range of values (much wider than what reported 
in the plots), in terms of the maximal converged Choi 
fidelity F, and to establish that (72 — 2.2 yields the clos- 
est results to an ideal CZ gate. Clearly, in practice, such 
couplings will not always be tunable at will, or possibly 
only within a given windows of values: it is anyway re- 
markable to be able to identify optimal values given a 
specific dynamical figure of merit. 

Let us now move on to a case which is not analytically 
treatable, with e = A = 1 and co m = 0.1m for 1 < m < 
10. As usual, we set g\ = 1 and consider different values 
of gi and different temperatures (zero temperature and 
f3 = 10): the values of F for such a configuration are 
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FIG. 9: MCE results for the Choi fidelity F versus rescaled 
time, for H rw with e = A = 1, gi = 1, M = 10, uj m = 0.1m 
for 1 < m < 10, /3 = 10 and different values of g2- The line 
F = 0.25 is reported for reference. 



displayed in Figs. [7]fT0l The initial peak in F is still 
apparent, but the breaking of the phase invariance by 
the term A clearly degrades the quality of the gate, with 
a maximum Choi fidelity which is now around 0.7, even 
at zero temperature. 

We now turn to the full spin-boson like Hamiltonian 
H (including the counter-rotating terms in the qubits- 
field coupling), and consider the case e = A = 1 and 
uj m = 0.1m for 1 < m < 10. The inclusion of the 
counter-rotating terms makes the simulation much more 
challenging to run and converge. Roughly speaking, the 
main difficulty one encounters comes down to the fact 
that the time-derivatives of the phase space positions of 
the basis grid, determined by the Ehrenfest dynamics as 
per Eq. ((6]), are much larger if the counter- rotating terms 



7 




ill ^^IHII 

IP 

III 1 IB 



II _ 

III B 
HlIRI I 



FIG. 10: MCE results for the Choi fidelity F versus rescaled 
time, for H rw with e = A = 1, g\ = 1, M = 10, cj m = 0.1m 
for 1 < m < 10, j3 = 10 and different values of g-z (red stands 
for higher values, blue for lower values). 
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FIG. 12: MCE results for the Choi fidelity F versus rescaled 
time, for H with e = A = 1, gi = 0.5, M = 10, cj m = 0.1m 
for 1 < m < 10, P = 10 and different values of $2- The line 
F = 0.25 is reported for reference. 




FIG. 11: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H with s = A = 1, gi = 0.5, 
M = 10, ui m = 0.1m for 1 < m < 10 and different values of 
Q2- The line F — 0.25 is reported for reference. 
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FIG. 13: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H rw with e = A = 1, gi = 1, 
M = 20, uim = 0.1m for 1 < m < 10, u m = 0.1(m - 10) for 
11 < m < 20 and different values of g2- The line F = 0.25 is 
reported for reference. 



are included. The time-dependent grid thus evolves much 
more rapidly in phase space and is likely to leave the 
dynamically relevant region and accumulate substantial 
errors earlier. We were however able to obtained well con- 
verged results by reducing the coupling g± to 0.5, which 
is still very far from the perturbativc regime. Likewise, 
we scanned values of gi up to 0.5. Quite significantly 
- as confirmed by Figs. [TT] and [121 respectively for zero 
temperature and /5 = 10 — we could not find any value 
of (72 such that the Choi fidelity of the CZ gate reaches 
0.25. In fact, strong enough couplings are necessary to 
entangle the two qubits on short enough time-scales but, 
with such strong couplings, the counter-rotating terms 
heat the qubits up too quickly for coherent effects to 
take place, at least in this region of parameters. This 
heating overshadows the effect of thermal fluctuations in 



the field, which are barely noticeable for j3 = 10 (and are 
instead manifest in the rotating wave regime at the same 
temperature) . 

In order to simulate the effect of a band of a 1- 
dimcnsional photonic band-gap medium, where modes 
are usually doubly degenerate in frequency (since they 
can propagate in either spatial direction), we have also 
considered a case with M — 20 modes, two for each 
equally spaced frequency. All the other parameters have 
been kept as above, with £ = A = 1 and g\ = 1 in the 
rotating wave case (Fig. IT3")) . and g\ = 0.5 in the full 
Hamiltonian (Fig. [H]). Comparing Fig. [13] with Fig. [7] 
shows that the initial peak in Choi fidelity is still present: 
moreover, not only does it occur earlier by a factor y2 (as 
expected because of the cooperation between the modes 
due to the balanced coupling), but it is also higher. Con- 
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FIG. 14: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H with e = A = 1, g\ = 0.5, 
M = 20, Lu m = 0.1m for 1 < m < 10, u) m = 0.1(m - 10) for 
11 < m < 20 and different values of g2- The line F = 0.25 is 
reported for reference. 



FIG. 16: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H rw with e — A — 1 and a com- 
mon Ohmic bath with a = 0.09, different cutoff frequencies 
and different numbers of total bath modes. The line F = 0.25 
is reported for reference. 
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FIG. 15: MCE results for the Choi fidelity F versus rescaled 
time at zero temperature, for H with e = A = 1 and a com- 
mon Ohmic bath with a — 0.09, uj c = 2.5 and different num- 
bers of total bath modes. The line F = 0.25 is reported for 
reference. 



trary to common intuition, this example shows that a 
larger number of mediating modes, in favourable dynam- 
ical configurations such as this, can actually be advan- 
tageous for the implementation of locally coherent dy- 
namics. Note that, for M = 20 modes, we needed about 
N = 400 coupled coherent states to achieve converged 
results. This is as large a basis set as we used in this 
study. 



A. Zero temperature Ohmic spin-boson bath 

The notion of entangling separated systems and of dis- 
tributing quantum coherence by interaction with com- 
mon heat baths or other incoherent means is well estab- 
lished in the quantum information and condensed matter 



communities, and has been explored under a number of 
- cither more specific and applied or more general and 



abstract - viewpoints [5 



However, the problem of 



studying the non-perturbative interaction of two qubits 
with a common bath is still in general a difficult one. 
Thus, to provide the reader with further evidence of the 
versatility and power of our approach, wc also report the 
application of the MCE method to the controlled-Z Choi 
fidelity for the case of the spin boson Hamiltonian H, 
with e = A = 1, and both qubits interacting with a com- 
mon bath at zero temperature and with Ohmic spectral 
density J(w) given by 



J(w) = —awe 

7T 



(22) 



where a is the Kondo parameter, which we fix at 0.09, 
and w c is a cutoff frequency. 

We use a standard approach to discretise the bath, 
which has already been proven very reliable for single 
spin Ohmic spin-boson systems [4lj . In particular, the 
frequencies and coupling strengths are chosen as follows: 



= — Lu r In 



771 ( 1 



1 - 



M 



(23) 



1 2 A 

9m 9m \ 



oJmauJc 1 — e 



2M 



(24) 



where U3 max is a free parameter of the numerics, which 
we will converge our results against (in that we choose it 
large enough to obtain converged results). In particular, 
we choose uj max = 12.5 for lo c = 2.5 and uj ma x = 6 for 
lu c = 1. Also, the coupling strengths g 3 m are defined as in 
Eqs. fllT J) and (fT2 |) . 
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FIG. 17: MCE results for the concurrence versus rescaled time at different temperatures for H rm with gi = 1 and gi = 2.1. In 
(a), e = A = 0, M = 10 (with w m = 0.1m for 1 < m < 10) and the initial state is | 4) = | f\); in (b), £ = A = 0, M = 10 
(with uj m — 0.1m for 1 < m < 10) and the initial state is | 2) = | i n ( c )i e = A = 1, the initial state is | 2) = | ft) and, 
respectively, M = 10 (with uj m — 0.1m for 1 < m < 10) for the dash dotted line and M — 20 (with ui m = 0.1m for 1 < m < 10 
and u: m = 0.1 (m — 10) for 11 < m < 20) for the dashed line. 



Fig. [TS] shows the convergence of our results for uj c = 
2.5 and the full Hamiltonian H in terms of the total num- 
ber of modes in the bath (M = 50 and M — 100). Quite 
interestingly, if the counter-rotating terms are included 
we obtain larger gate fidelities when mimicking a bath 
than for a smaller set of discrete bus frequencies. Off 
resonant modes are more influential in the full Hamilto- 
nian and seems to be captured faithfully by our method. 

In Fig. [THl instead, we report results for the rotating 
wave Hamiltonian H rw and different cutoff frequencies 
and number of bath modes. In this case, the effect of 
the counter-rotating terms is rather limited. Also, the 
influence of larger cutoff frequencies at longer times is 
clearly visible in the plot. 



VI. CONCLUSIONS 



We have presented an extensive numerical study, based 
on Multi-Configurational Ehrenfest trajectories, on the 
dynamics of two qubits interacting with a common set of 
bosonic field modes, obtaining converged results for the 
Choi fidelity of an entangling CZ gate between the qubits 
for a rather wide range of Hamiltonian parameters and 
field temperatures, which cannot be covered by pertur- 
bation theory or other approximate approaches. We thus 
demonstrated the capability of tracking, analysing in de- 
tail, and even optimising with respect of certain ranges of 
some parameters, specific aspects of the coherent quan- 
tum dynamics of the qubits. 



V. ENTANGLEMENT GENERATION 

Typically, a large Choi fidelity for the (entangling) CZ 
gate corresponds to the generation of substantial entan- 
glement between the two qubits. To support this state- 
ment, we report here a brief study on the entanglement 
generated between two qubits. As an entanglement quan- 
tifier, we adopt the concurrence, an entanglement mono- 
tone that can be easily calculated for a system of two 
qubits [67|. Figs. 17(a)| 17(b) and 17(c) show the concur- 
rence versus rescaled time for the rotating wave Hamil- 
tonian H rw with different initial states, temperatures, 
dynamical parameters and number of modes. The degra- 
dation of quantum entanglement due to temperature is 



apparent (Figs 17(a) and 17(b) ), along with the speed up 



in the entanglement generation induced by a doubling of 
the modes (Fig. 17(c)). 

It is also worth noticing that we did not find any 
region of parameters where entanglement between the 
two qubits is generated for the full Hamiltonian H and 
AI = 10, thus mirroring our failure in obtaining Cz fi- 
delities larger than 0.25. 



We were able to properly take into account the effect 
of finite bath's temperatures on the reduced dynamics 
of the qubits, and also to highlight some counterintuitive 
features related to the scaling of coherent signatures with 
the number of field modes (which we varied over the 
range 1 — 100), showing that at times more mediating 
modes can actually be advantageous for the distribution 
of quantum coherence. 

The main limitations of our approach lie in the dif- 
ficulty of handling counter-rotating qubit-field coupling 
terms in the strong coupling regime (i.e., when the cou- 
pling strengths are comparable to the inherent dynami- 
cal frequencies of the qubits). Even in such instances, we 
could however reach convergence by somewhat limiting 
the range of the coupling strengths. 

Within such limitations, the MCE approach has hence 
been established as a powerful tool for the detailed study 
of complex quantum dynamics even with relatively lim- 
ited resources (desktop computers), typically for systems 
where discrete sets of up to 100 field modes are involved. 
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FIG. 18: Norm and expectation value of the energy for MCE 
results at zero temperature, for H rw with e = A = gi = 1, 
g 2 = 2.7, M = 10 (with u m = 0.1m for 1 < m < 10) and 
different values of TV and comp. 
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FIG. 21: Entries of the qubits' density matrix gn and Q13 
for MCE results at zero temperature, for H with e — A = 
gi = g2 = 1, M = 10 (with LJ m = 0.1m for 1 < m < 10) and 
different values of N and comp. 
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Appendix A: Ehrenfest dynamics of the coherent 
states 



Here, we elaborate on Eq. ^ and derive explicitly the 

equation of motion of each complex parameter a^" 1 ' • The 
approximated Lagrangian Cj reads 



FIG. 19: Entries of the qubits' density matrix gn and gi3 for 
MCE results at zero temperature, for H rw with e = A = g\ = 
1, 32 = 2.7, M = 10 (with u m = 0.1m for 1 < m < 10) and 
different values of N and comp. 
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FIG. 20: Norm and expectation value of the energy for MCE 
results at zero temperature, for H with e = A = gi = g% = 1, 
M = 10 (with uj m = 0.1m for 1 < m < 10) and different 
values of N and comp. 



Cj = c *,j c n,j{<XjJ\H\<Xj,n) -»X) c *.Aj 

l,n=l 1 = 1 

-it\^e(^.-^), (ad 



1=1 



where 1 1, otj) = \ I) ® | ctj). Note that (ctj,l \H\ atj,n) is 
just a function of the vector ctj and the integers I and n, 
promptly evaluated by normal ordering H . 
The Euler-Lagrange equation for a*:" 1 ' then is: 



da 



(m) 



i,n=l 9a j 

d • (m)* 

+*Ei^ 12 j 



1=1 



dt 



2 

Em 



(A2) 



d (JL a {m) * 



d__ac z _ 
d7^~p 
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which, by neglecting the time dependence of c/j (setting 



d_ ( 

dt \l^i=i 



=1 Q, 



0) , can be rearranged to obtain 



a) =1 



r n <3 Sa (m) 



[oij, I \H\ ttj, n) 



Ei=i l c ul 



(A3) 



This equation governs the evolution of the vectors OLj in 
our numerics. 
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FIG. 22: Choi fidelity and concurrence for two different separable initial states versus rescaled time, for H rm with e — A = 0, 
gi = 1, g<x = 2.1, ft = 0.5 and different values of Nt and M. In all plots M = 10 with uj m = 0.1m for 1 < m < 10. "Conjugate" 
refers to the fact that for those curves the initial centres of the coherent states ctj are in complex conjugate pairs. 



Appendix B: Dynamics of the state vector 

For the sake of completeness, let us also report the sys- 
tem of dynamical equations (Schrodinger equation on the 
subspace spanned by the basis grid) for the parameters 
qj, which can be derived by Eq. (|2|) and reads 



N 

E 

fc=l L 



/ x a j a j a j a j 
i^jkCLk + iQjk I ajdfe J - J -^— 

d 

-y](aj,l \H\ h, a. k )VL. jk c h ,k 
h=i 



where = {oLj |atfc). To enhance the stability of the 
numerical treatment, the parameters qj are actually re- 
defined by multiplication with a smooth phase factor (es- 
sentially a semi-classical action). 



Appendix C: Convergence of MCE results 

To give an idea of the quality and range of reliabil- 
ity of our results, we provide here some evidence of the 
convergence of our numerics. 

Throughout our study, the centres of the initial set 
of coherent states are distributed in phase space with a 
Gaussian distribution with standard deviation l/comp. 
The parameter camp is a free parameter of the study, 



which is tuned to optimise convergence. As indicators 
of the quality of the numerics we will observe the con- 
vergence of specific entries of the density matrix of the 
two qubits g, as well as the 'norm' Tr(p) and the expec- 
tation value of the energy Tr(iJg), which are obviously 
conserved in the exact dynamics. Notice that our method 
does not have any in-built routine guaranteeing the con- 
servation of the state vector's norm, so that Ti(g) is a 
relevant figure of merit to assess its reliability. 



Figs. 18(a) and 18(b) display the norm and expectation 



value of the energy for a case of non-number conserving 
rotating wave Hamiltonian, while in Figs. 19(a) and |19(bj] 
the entries Qn and Q\z are plotted. The reliability of the 
numerics over the whole timeframe considered is appar- 
ent (for large enough compression parameter comp), in 
terms of both convergence with increasing number ./V of 
coherent states and of conservation of invariant quanti- 
ties. As anticipated, the situation is much more dire for 
the fu ll H amilton ian H. In this case, Figs. |20(a)| |20(b)| 
21(a) and 21(b)) show that our numerics are only reli- 



able up to rescaled times around 2.5, after which both 
convergence, and norm and energy conservation are lost, 
even at smaller coupling strengths (in that <?2 = 1 rather 
than 52 = 2.7 as before). Finally, we show three exam- 
ples of convergence of our results at finite temperature 
(here, /3 = 0.5) with respect to the increase in the num- 
ber of states Nt over which the thermal distribution of 
Eq. HE} is sampled (Figs. |22(a)||22(b)| and |22(c) |. 
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